Lorenzo Biasi and Michael Aichmüller
In [2]:
import numpy as np
from scipy.special import binom
import matplotlib.pylab as plt
from scipy.misc import factorial as fact
%matplotlib inline
def binomial(p, n, k):
return binom(n, k) * p ** k * (1 - p) ** (n-k)
$\Omega$ will be all the possible combinations we have for 150 object two have two diffent values. For example (0, 0, ..., 0), (1, 0, ..., 0), (0, 1, ..., 0), ... (1, 1, ..., 0), ... (1, 1, ..., 1). This sample space has size of $2^{150}$. The random variable $X(\omega)$ will be the number of defective objects there are in the sample $\omega$. We can also define $Y(\omega) = 150 - X(\omega)$, that will be counting the number of checked items.
In [5]:
p = 1. / 365
1 - np.sum(binomial(p, 23 * (22) / 2, 0))
Out[5]:
In [3]:
np.sum(binomial(p, 150, np.arange(5, 10)))
Out[3]:
In [4]:
plt.bar(np.arange(20), binomial(p, 150, np.arange(20)))
plt.bar(np.arange(5), binomial(p, 150, np.arange(5)))
plt.bar(np.arange(5, 10), binomial(p, 150, np.arange(5,10)))
plt.xlabel('# defectives')
plt.ylabel('P(X=k)')
Out[4]:
For computing how big $q$ needs to be we can compute the probability $p^*$ that nobody has the same birthday in a group of $q$ and compute $1 - p^*$. The first two people will not have the same birthday with probability of $364/365$, the probability that the third will also have a different birthday will be $364/365 * 363 / 365$. this will go on until the last person. One can make the computation and finds that the minimum for having over 50% of probability that at least two people have the same birthday is 23 with p = 50.73%.
In [5]:
def not_same_birthday(q):
return np.prod((365 - np.arange(q))/ 365)
q = 45
p = np.empty(q - 1)
for i in range(1, q):
p[i - 1] = 1 - not_same_birthday(i)
plt.plot(np.arange(1, q), p)
plt.plot(23, 1 - not_same_birthday(23), 'r+', label='23 people')
plt.grid()
plt.ylabel('Probability')
plt.xlabel('q')
plt.legend()
1 - not_same_birthday(23)
Out[5]:
In [6]:
import itertools
x = [1, 2, 3, 4, 5, 6]
omega = set([p for p in itertools.product(x, repeat=3)])
print(r'Omega has', len(omega), 'elements and they are:')
print(omega)
X would be -30 when the sample $\omega$ has no 6s, 50 when has one, 75 when it has two, and 100 when it has three. The probability distribution of such variable would be the binomial with $p = 1 / 6$, $n=3$ and $k$ the number of 6s.
So:
$P_X(X = -30) = {3\choose 0}(1 / 6)^0(1-1/6)^{3-0}$
$P_X(X = 50) = {3\choose 1}(1 / 6)^1(1-1/6)^{3-1}$
$P_X(X = 75) = {3\choose 2}(1 / 6)^2(1-1/6)^{3-2}$
$P_X(X = 100) = {3\choose 3}(1 / 6)^3(1-1/6)^{3-3}$
In [7]:
g = binomial(1 / 6, 3, np.arange(4)) * np.array([-30, 50, 75, 100])
np.sum(g)
Out[7]:
In [8]:
plt.bar(np.arange(4), g)
plt.plot([-.5, 3.5], np.ones(2) * np.sum(g), 'r')
Out[8]:
In [ ]: